{
 "metadata": {
  "signature": "sha256:2baa25856855a6270066c29406122ce9263cf1916bab38630886d2a53a5f499e"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Matplotlib: gridding irregularly spaced data\n",
      "======================================================================\n",
      "\n",
      "A commonly asked question on the matplotlib mailing lists is \"how do I\n",
      "make a contour plot of my irregularly spaced data?\". The answer is,\n",
      "first you interpolate it to a regular grid. As of version 0.98.3,\n",
      "matplotlib provides a griddata function that behaves similarly to the\n",
      "matlab version. It performs \"natural neighbor interpolation\" of\n",
      "irregularly spaced data a regular grid, which you can then plot with\n",
      "contour, imshow or pcolor.\n",
      "\n",
      "Example 1\n",
      "---------\n",
      "\n",
      "This requires Scipy 0.9:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from scipy.interpolate import griddata\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy.ma as ma\n",
      "from numpy.random import uniform, seed\n",
      "# make up some randomly distributed data\n",
      "seed(1234)\n",
      "npts = 200\n",
      "x = uniform(-2,2,npts)\n",
      "y = uniform(-2,2,npts)\n",
      "z = x*np.exp(-x**2-y**2)\n",
      "# define grid.\n",
      "xi = np.linspace(-2.1,2.1,100)\n",
      "yi = np.linspace(-2.1,2.1,100)\n",
      "# grid the data.\n",
      "zi = griddata((x, y), z, (xi[None,:], yi[:,None]), method='cubic')\n",
      "# contour the gridded data, plotting dots at the randomly spaced data points.\n",
      "CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')\n",
      "CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.jet)\n",
      "plt.colorbar() # draw colorbar\n",
      "# plot data points.\n",
      "plt.scatter(x,y,marker='o',c='b',s=5)\n",
      "plt.xlim(-2,2)\n",
      "plt.ylim(-2,2)\n",
      "plt.title('griddata test (%d points)' % npts)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "![](files/attachments/Matplotlib_Gridding_irregularly_spaced_data/griddataexample1.png)\n",
      "\n",
      "Example 2\n",
      "---------"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from matplotlib.mlab import griddata\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy.ma as ma\n",
      "from numpy.random import uniform\n",
      "# make up some randomly distributed data\n",
      "npts = 200\n",
      "x = uniform(-2,2,npts)\n",
      "y = uniform(-2,2,npts)\n",
      "z = x*np.exp(-x**2-y**2)\n",
      "# define grid.\n",
      "xi = np.linspace(-2.1,2.1,100)\n",
      "yi = np.linspace(-2.1,2.1,100)\n",
      "# grid the data.\n",
      "zi = griddata(x,y,z,xi,yi)\n",
      "# contour the gridded data, plotting dots at the randomly spaced data points.\n",
      "CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')\n",
      "CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.jet)\n",
      "plt.colorbar() # draw colorbar\n",
      "# plot data points.\n",
      "plt.scatter(x,y,marker='o',c='b',s=5)\n",
      "plt.xlim(-2,2)\n",
      "plt.ylim(-2,2)\n",
      "plt.title('griddata test (%d points)' % npts)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "By default, griddata uses the scikits delaunay package (included in\n",
      "matplotlib) to do the natural neighbor interpolation. Unfortunately, the\n",
      "delaunay package is known to fail for some nearly pathological cases. If\n",
      "you run into one of those cases, you can install the matplotlib natgrid\n",
      "toolkit. Once that is installed, the griddata function will use it\n",
      "instead of delaunay to do the interpolation. The natgrid algorithm is a\n",
      "bit more robust, but cannot be included in matplotlib proper because of\n",
      "licensing issues.\n",
      "\n",
      "The radial basis function module in the scipy sandbox can also be used\n",
      "to interpolate/smooth scattered data in n dimensions. See\n",
      "[\"Cookbook/RadialBasisFunctions\"] for details.\n",
      "\n",
      "Example 3\n",
      "---------\n",
      "\n",
      "A less robust but perhaps more intuitive method is presented in the code\n",
      "below. This function takes three 1D arrays, namely two independent data\n",
      "arrays and one dependent data array and bins them into a 2D grid. On top\n",
      "of that, the code also returns two other grids, one where each binned\n",
      "value represents the number of points in that bin and another in which\n",
      "each bin contains the indices of the original dependent array which are\n",
      "contained in that bin. These can be further used for interpolation\n",
      "between bins if necessary.\n",
      "\n",
      "The is essentially an Occam's Razor approach to the matplotlib.mlab\n",
      "griddata function, as both produce similar results."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# griddata.py - 2010-07-11 ccampo\n",
      "import numpy as np\n",
      "\n",
      "def griddata(x, y, z, binsize=0.01, retbin=True, retloc=True):\n",
      "    \"\"\"\n",
      "    Place unevenly spaced 2D data on a grid by 2D binning (nearest\n",
      "    neighbor interpolation).\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    x : ndarray (1D)\n",
      "        The idependent data x-axis of the grid.\n",
      "    y : ndarray (1D)\n",
      "        The idependent data y-axis of the grid.\n",
      "    z : ndarray (1D)\n",
      "        The dependent data in the form z = f(x,y).\n",
      "    binsize : scalar, optional\n",
      "        The full width and height of each bin on the grid.  If each\n",
      "        bin is a cube, then this is the x and y dimension.  This is\n",
      "        the step in both directions, x and y. Defaults to 0.01.\n",
      "    retbin : boolean, optional\n",
      "        Function returns `bins` variable (see below for description)\n",
      "        if set to True.  Defaults to True.\n",
      "    retloc : boolean, optional\n",
      "        Function returns `wherebins` variable (see below for description)\n",
      "        if set to True.  Defaults to True.\n",
      "   \n",
      "    Returns\n",
      "    -------\n",
      "    grid : ndarray (2D)\n",
      "        The evenly gridded data.  The value of each cell is the median\n",
      "        value of the contents of the bin.\n",
      "    bins : ndarray (2D)\n",
      "        A grid the same shape as `grid`, except the value of each cell\n",
      "        is the number of points in that bin.  Returns only if\n",
      "        `retbin` is set to True.\n",
      "    wherebin : list (2D)\n",
      "        A 2D list the same shape as `grid` and `bins` where each cell\n",
      "        contains the indicies of `z` which contain the values stored\n",
      "        in the particular bin.\n",
      "\n",
      "    Revisions\n",
      "    ---------\n",
      "    2010-07-11  ccampo  Initial version\n",
      "    \"\"\"\n",
      "    # get extrema values.\n",
      "    xmin, xmax = x.min(), x.max()\n",
      "    ymin, ymax = y.min(), y.max()\n",
      "\n",
      "    # make coordinate arrays.\n",
      "    xi      = np.arange(xmin, xmax+binsize, binsize)\n",
      "    yi      = np.arange(ymin, ymax+binsize, binsize)\n",
      "    xi, yi = np.meshgrid(xi,yi)\n",
      "    \n",
      "    # make the grid.\n",
      "    grid           = np.zeros(xi.shape, dtype=x.dtype)\n",
      "    nrow, ncol = grid.shape\n",
      "    if retbin: bins = np.copy(grid)\n",
      "    \n",
      "    # create list in same shape as grid to store indices\n",
      "    if retloc:\n",
      "        wherebin = np.copy(grid)\n",
      "        wherebin = wherebin.tolist()\n",
      "\n",
      "    # fill in the grid.\n",
      "    for row in range(nrow):\n",
      "        for col in range(ncol):\n",
      "            xc = xi[row, col]    # x coordinate.\n",
      "            yc = yi[row, col]    # y coordinate.\n",
      "    \n",
      "            # find the position that xc and yc correspond to.\n",
      "            posx = np.abs(x - xc)\n",
      "            posy = np.abs(y - yc)\n",
      "            ibin = np.logical_and(posx < binsize/2., posy < binsize/2.)\n",
      "            ind  = np.where(ibin == True)[0]\n",
      "            \n",
      "            # fill the bin.\n",
      "            bin = z[ibin]\n",
      "            if retloc: wherebin[row][col] = ind\n",
      "            if retbin: bins[row, col] = bin.size\n",
      "            if bin.size != 0:\n",
      "                binval         = np.median(bin)\n",
      "                grid[row, col] = binval\n",
      "            else:\n",
      "                grid[row, col] = np.nan   # fill empty bins with nans.\n",
      "    \n",
      "    # return the grid\n",
      "    if retbin:\n",
      "        if retloc:\n",
      "            return grid, bins, wherebin\n",
      "        else:\n",
      "            return grid, bins\n",
      "    else:\n",
      "        if retloc:\n",
      "            return grid, wherebin\n",
      "        else:\n",
      "            return grid"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The following example demonstrates a usage of this method."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import griddata\n",
      "\n",
      "npr = np.random\n",
      "npts = 3000.                            # the total number of data points.\n",
      "x = npr.normal(size=npts)            # create some normally distributed dependent data in x.\n",
      "y = npr.normal(size=npts)            # ... do the same for y.\n",
      "zorig = x**2 + y**2                      # z is a function of the form z = f(x, y).\n",
      "noise = npr.normal(scale=1.0, size=npts) # add a good amount of noise\n",
      "z = zorig + noise                    # z = f(x, y) = x**2 + y**2\n",
      "\n",
      "# plot some profiles / cross-sections for some visualization.  our\n",
      "# function is a symmetric, upward opening paraboloid z = x**2 + y**2.\n",
      "# We expect it to be symmetric about and and y, attain a minimum on\n",
      "# the origin and display minor Gaussian noise.\n",
      "\n",
      "plt.ion()   # pyplot interactive mode on\n",
      "\n",
      "# x vs z cross-section.  notice the noise.\n",
      "plt.plot(x, z, '.')\n",
      "plt.title('X vs Z=F(X,Y=constant)')\n",
      "plt.xlabel('X')\n",
      "plt.ylabel('Z')\n",
      "\n",
      "# y vs z cross-section.  notice the noise.\n",
      "plt.plot(y, z, '.')\n",
      "plt.title('Y vs Z=F(Y,X=constant)')\n",
      "plt.xlabel('Y')\n",
      "plt.ylabel('Z')\n",
      "\n",
      "# now show the dependent data (x vs y).  we could represent the z data\n",
      "# as a third axis by either a 3d plot or contour plot, but we need to\n",
      "# grid it first....\n",
      "plt.plot(x, y, '.')\n",
      "plt.title('X vs Y')\n",
      "plt.xlabel('X')\n",
      "plt.ylabel('Y')\n",
      "\n",
      "# enter the gridding.  imagine drawing a symmetrical grid over the\n",
      "# plot above.  the binsize is the width and height of one of the grid\n",
      "# cells, or bins in units of x and y.\n",
      "binsize = 0.3\n",
      "grid, bins, binloc = griddata.griddata(x, y, z, binsize=binsize)  # see this routine's docstring\n",
      "\n",
      "\n",
      "# minimum values for colorbar. filter our nans which are in the grid\n",
      "zmin    = grid[np.where(np.isnan(grid) == False)].min()\n",
      "zmax    = grid[np.where(np.isnan(grid) == False)].max()\n",
      "\n",
      "# colorbar stuff\n",
      "palette = plt.matplotlib.colors.LinearSegmentedColormap('jet3',plt.cm.datad['jet'],2048)\n",
      "palette.set_under(alpha=0.0)\n",
      "\n",
      "# plot the results.  first plot is x, y vs z, where z is a filled level plot.\n",
      "extent = (x.min(), x.max(), y.min(), y.max()) # extent of the plot\n",
      "plt.subplot(1, 2, 1)\n",
      "plt.imshow(grid, extent=extent, cmap=palette, origin='lower', vmin=zmin, vmax=zmax, aspect='auto', interpolation='bilinear')\n",
      "plt.xlabel('X values')\n",
      "plt.ylabel('Y values')\n",
      "plt.title('Z = F(X, Y)')\n",
      "plt.colorbar()\n",
      "\n",
      "# now show the number of points in each bin.  since the independent data are\n",
      "# Gaussian distributed, we expect a 2D Gaussian.\n",
      "plt.subplot(1, 2, 2)\n",
      "plt.imshow(bins, extent=extent, cmap=palette, origin='lower', vmin=0, vmax=bins.max(), aspect='auto', interpolation='bilinear')\n",
      "plt.xlabel('X values')\n",
      "plt.ylabel('Y values')\n",
      "plt.title('X, Y vs The No. of Pts Per Bin')\n",
      "plt.colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The binned data:\n",
      "![](files/attachments/Matplotlib_Gridding_irregularly_spaced_data/bin_small.png)\n",
      "\n",
      "Raw data superimposed on top of binned data:\n",
      "![](files/attachments/Matplotlib_Gridding_irregularly_spaced_data/raw_small.png)"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}